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ABSTRACT 



The use of the Sequential Regression Algorithm (SER) to 
coherently remove background noise from an ELF sensor is presented. 
The SER algorithm is described for a multi-channel application in 
order to cancel coherent portions of reference sensors from a 
primary sensor. The algorithm adaptively accounts for differences 
between two parallel array platforms for the purpose of coherent 
subtraction. A section on likelihood detector schemes is also 
presented. This work is in support of a submerged ELF sensor array 
project run by the Johns Hopkins University Applied Physics Lab. 
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I . INTRODUCTION 

A. BACKGROUND 

Extremely low frequency (ELF) geoelectric noise has a 
highly impulsive nature which makes it very difficult to 
remove by filtering. The main contributors to the 
impulsiveness of the ELF spectrum are lightning discharges. 
The lowest frequencies of the ELF range (1 to 60Hz) have such 
a low attenuation rate as they propagate, that lighting 
discharges all over the world must be considered when studying 
the noise spectrum at a single location. ELF signals have 
been studied for such applications as power transfer, 
communications, geophysical surveys, and so forth. The amount 
of information on the subject is immense and increasing. How 
to account for the impulsive background electromagnetic fields 
inevitably present in any system operating in this frequency 
range remains a serious problem that any system designer must 
answer . 

The background noise of the lower ELF frequency range can 
be thought of as consisting of both local and global 
phenomena. The global component is comprised of signals that 
are relatively coherent over extended distances. A detailed 
recording of data in this frequency range 

...is characterized by bursts of a few cycles duration 
which may be ten or more times the mean amplitude between 
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bursts. This bursty character is apparently coherent over 
most, if not all, of the Earth. [Ref. 1] 

The local phenomenon is comprised of those components that are 

not coherent over extended distances and are assumed to be 

products of the local environment of the detector. A set of 

parallel detectors separated by some distance can be 

coherently combined to remove the global portion of the 

background noise and thus lower the noise floor. 

B . OBJECTIVE 

The Applied Physics Lab (APL) at Johns Hopkins University 
has for years studied and made measurements of the background 
geoelectric noise on widely separated pairs of electrodes. 
The purpose of this thesis is to describe techniques of 
adaptive filtering that the author applied to data collected 
by APL, in order to determine the amount of noise 
cancellation. The techniques are implemented using software 
developed during an extended visit to the Applied Physics Lab 
by the author. In a related effort, software will also be 
presented implementing appropriate detection methodology for 
the detection of narrowband signals embedded in the noise. 

C . ORGANIZATION 

The thesis begins with a brief discussion of the 
constituents of this ELF noise to demonstrate the difficulty 
in modeling of the spectrum and to provide supporting 
background information. The basic concepts of the adaptive 
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filtering will be presented as a lead in to the adaptive 
algorithm to be used to cancel the noise. These introductory 
sections will be followed by a discussion of stability and 
convergence concerns for the adaptive algorithm in this type 
of data environment. The detection discussion is then 
presented as a continuation of the processing required for the 
detection of narrowband signals. The adaptive algorithm is 
then used in the detector scheme to demonstrate performance of 
the developed software. A final section will summarize the 
work performed, draw relevant conclusions about the work and 
propose possible future studies that could be derived from 
this work. 
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II. THE ELF BACKGROUND 



Many studies have been conducted in an attempt to develop 
an understanding of what actually constitutes the natural and 
man-made background noise spectrum in the one to sixty hertz 
range. The sources of the energy in this frequency range can 
be primarily attributed to atmospheric disturbances and man- 
made power transfer systems. A major problem with these 
sources of noise is that they are highly irregular in both 
magnitude and duration of the noise signal. Attempts to 
classify and assign specific characteristics to these sources 
have proven difficult. A brief history of the research 
conducted in these areas is provided in the following two 
subsections . 

A. ATMOSPHERIC NOISE 

The source of the majority of noise in the one to sixty 
hertz frequency range is lightning discharges. The lower ELF 
range discussed here also contains what is called the Schumann 
Range. The Schumann range takes its name from the early 
theoretical work of the German scientist W. 0. Schumann. 
Schumann theorized that the cavity created by the surface of 
the Earth and the Ionosphere has naturally occurring resonant 
frequencies. He took this research further to derive a set of 
harmonic resonant equations for calculating these frequencies. 
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The work of Schumann has been carried further with the 



actual measurement of the Schumann resonance frequencies in 
field experiments. In 1959 and 1960, H. L. von Konig 
substantiated the presence of Schumann resonances by observing 
the waveforms in the output of a narrow band amplifier. The 
experimental presence of Schumann resonances created a large 
quantity of literature on the subject as scientists attempted 
to understand the phenomenon. Once a foundation of the 
propagation properties of the Earth-Ionosphere cavity was 
established, the search for the source of the excitation 
began . [Ref . 2 ] 

The most prominent excitation source of the Earth- 
Ionosphere cavity is the cloud to ground lightning stroke. 
The number of similar strokes present in any one flash 
observed during a thunderstorm is a random variable that 
normally ranges from two to twelve. Since the number of 
strokes is generally too high to measure individual waveforms, 
a more applicable measure is the power spectral 
density [Ref. 3]. The frequency spectrum of this noise was 
resolved in sufficient detail to estimate the quality factor, 
Q, of the earth-ionosphere cavity by M. Balser and C. A. 
Wagner in 1960. The Q of a cavity resonator is a measure of 
the bandwidth of the resonator [Ref . 4] . The source of the 
noise, as theorized by H. Raemer in 1961, is the response of 
the earth-ionosphere cavity to electrical discharges created 
in thunderstorms all over the world. Raemer attempted to 



5 



model the ionosphere in an attempt to reproduce mathematically 
the measured response of Balser and Wagner. Although Raemer ' s 
attempt fell short of its goal, it began the process of 
refinement of a working model of the ionosphere that continues 
today . [Ref . 2 ] 

The early work performed by these scientists brought out 
the difficulty in attempting to model this naturally occurring 
phenomenon due to its randomness. Parameters must be 
estimated by some method in order to model this random 
process. The highly variable nature of these parameters makes 
the task of modeling the noise nearly impossible. Examples of 
these varying parameters are diurnal variations in the 
ionosphere, twenty-four hour variations in the ionosphere, 
seasonal variations in the locations of thunderstorm regions, 
and the randomness of individual lightning strikes. These are 
only a few examples of the parameters that must be considered 
in order to model the atmospheric noise. 

B. MAN-MADE NOISE 

The presence of atmospheric noise is not the only concern 
that must be addressed when discussing background noise in the 
one to sixty hertz range. Many components in today's world 
emit unshielded interference that falls in this range. These 
components are the result of either faulty equipment or poor 
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design. Examples of this noise source are power lines that 
are used to distribute alternating current through populated 
areas . 

Large ferric objects can create fields around them that 
can be detected. The strength of these fields can be well 
below the previously mentioned sources of noise but can be 
discernable in certain situations. Motion of these large 
magnetic objects in the Earth's natural magnetic field is the 
source of a small field. 

Motion of the detector can create an increase in the 
background noise. Any slight motion of the detector can 
create misleading information and in most cases, an increase 
in the noise signal. The background ELF spectrum can be 
thought of as the summation of three orthoganal vectors. 
Relative motion between these three vectors and the detector 
creates a false signal that can raise the noise floor. This 
type of noise can be compensated for by the use of motion 
sensor. [Ref . 5] 

Other sources of man-made noise can be found in the 
detector equipment itself. Improper shielding of cables and 
electronic devices used in the design of the detector and its 
processing components can lead to elevated noise floors. 
Unlike naturally occurring atmospheric noise, these sources 
can be traced and eliminated. The removal of this noise can 
come in two forms, elimination of the source or isolation of 
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the source from the detector. The fact that man-made noise 



sources exist must not be overlooked when attempting to 
process signals at the output of the detector. 
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III. ADAPTIVE FILTERS 



A. GENERAL TERMINOLOGY 

The ELF background noise has been described as a highly 
impulsive, non-stationary random process. Some of the causes 
of these variations have been briefly described above. The 
use of traditional filter designs may not be the optimum 
method for reduction of the noise floor for reception of low 
strength signals in the one to fifty hertz range. Since their 
introduction, the use of adaptive filters has shown promise in 
combating exactly this type of scenario. To understand the 
application of adaptive filters in this case, a brief 
discussion of relevant terms must be conducted. 

An adaptive filter is by definition a filter that has the 
capability to adjust, by itself, a set of design parameters 
that are based on estimated statistical characteristics of 
the signal to be filtered. This process of adjusting the 
filter to the present situation alleviates the need to have a 
priori knowledge of the relevant signal characteristics. The 
previous chapter described how difficult it is to derive a set 
of statistical characteristics for the ELF noise environment. 
If a priori knowledge of the relevant signal characteristics 
were known, an optimum filter could be designed. The goal of 
adaptive filtering is to estimate these statistical parameters 
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and to refine the estimation through the use of an algorithm. 
If an appropriate algorithm can be found, then after a 
sufficient number of iterations, the adaptive filter should 
converge to the optimum filter. The use of adaptive filters 
requires the designer to look for the optimum algorithm for 
estimating the relevant signal characteristics instead of 
attempting to discover the actual signal 
characteristics . [Ref • 6] 

Adaptive filters are divided into two basic categories, 
open-loop and closed-loop. An open-loop configuration is 
usually a two stage process. The first stage is used to 
"learn" the statistics of the relevant signal. The results of 
the first stage are then used in the second stage to compute 
the filter parameters required using a nonrecursive algorithm. 
In contrast, the closed-loop configuration uses only one stage 
to develop the filter parameters: the relevant signal 
characteristics are not explicitly estimated but are inferred 
from a recursive algorithm that updates the filter parameters 
directly as each new data point is obtained. The adaptive 
filter gains additional knowledge of the signal 
characteristics during each iteration. The gain in knowledge 
results in an improvement of the filter parameters which are 
adjusted and their performance is used as an input to the next 
iteration of the algorithm. The adaptive filter discussed in 
this thesis falls in the closed loop category. The advantage 
of closed-loop over open-loop is that the requirement for only 
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one stage generally corresponds to a less expensive 
configuration based on its simpler implementation. [Ref 6] 

The description of the filters that will be analyzed in 
this study carries terminology that is pertinent to all 
filters, not just adaptive filters. Filters are implemented 
in either continuous-time or discrete-time which defines the 
form of the relationship between the input and output of the 
filter. The filters that will be studied are of the discrete- 
time version which means that the filter may be described by 
a difference equation. The structure of the filters to be 
studied are referred to as finite-impulse response (FIR) , 
tapped-delay-line, or transversal filters. This description 
means that the filter's output relies only on the past and 
present values of the input. The advantage of a FIR filter is 
its inherent stability . [Ref 6] 

The algorithm that is used to update the filter parameters 



is generally 


named after some 


of its 


operating 


characteristics . 


The algorithm used 


here is 


called 


the 


Sequential Regression Algorithm or 


SER [Ref . 


7] . 


This 



algorithm uses an iterative approach to approximating the 
weighting factors required for a linear combiner. Figure 1 
shows a basic diagram of a linear combiner. 

This system uses simultaneous samples of the input signal 
at time index k, to produce an output signal, y k , by a linear 
operation with weighting factors, w 1 . A second way to 
physically interpret the linear combiner is to look at it in 
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the form of a transversal filter. Figure 2 shows how this 
system is designed: X k corresponds to the last L samples of 
x at time index k. The weighting factors are dual indexed on 




the time delay, 1, followed by the time index, k. The weight 
factor, w 2k , corresponds to the weight factor for the sample 
taken two samples ago from the present sample time, k. 
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B. THEORY AND DESIGN 



1. Configuration 

The filter to be developed here is a combination of 
the previous two figures. The transversal filter is applied 
to a set of multiple sensor outputs to form an input vector 
that is N x L data points long, where L is the number of 
separate inputs or sources, and N is the length of the 
transversal filter. There are two reasons for using a filter 
design of this nature. The first is the ability to use more 
than one source of input, allowing for a better spectral 
sampling of the background ELF energy. The second is the 
ability to use the previous N samples thereby smoothing the 
output since it is not dependent on instantaneous values which 
may vary greatly from sample to sample. Both of these ideas 
are beneficial to the overall performance of the filtering 
process . 

The output of the filter, y k , can be expressed in 
terms of the input matrix and the weighting factors matrix. 
To do this, some defining equation must be expressed for the 
input matrix, X k , 

x* - [* lt x 2k ... X u ]* (1) 

and weight factors matrix, W k , 

W 2 Jc ••• T < 2 > 

where each sensor's input vector is 
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X lk = [x^k) x 1 {k-l) x 1 {k-2) ... x 1 (k~(N-l))] 



( 3 ) 



and weighting factor vector is 

^lJt = ^KJc-l) *Uk-2) • • • W 1 (Jc-(N-l) ) 1 

The first subscript of Equations 3 and 4, labeled 1, 2, 3, 

. .., L, indicates the source of the sample. The second 
subscript, k, indicates the time of the sample. A total of N 
sample values from each sensor are used for each calculation 
of y k . A physical model for the filter to be studied is shown 
in Figure 3. The scalar output value, y k in Figure 3, is 
defined as: 

y k =xlw k =V T k X k . ( 5 ) 

These basic equations will become the starting point 
for the development of an operating filter system for the 
removal of the ELF background noise. It J s preferable that 
the weighting factors, W k , cause the output, y k , to be an 
estimate of some desired signal, d k . The concept of the 
filter is to use the weighting factors to model the 
differences between signals originating from the same source 
but sensed at various locations. The inputs, x 1 in Figure 3, 
represent samples of a signal taken at different locations. 
The desired signal, d k , is also a sample of the same 
background source signal. The difference between the desired 
signal and the output of the filter is expressed as the error 
in the estimate, e k . 
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*k = d k - y k 



( 6 ) 



In the system designed here, the estimate error is the 
output sought for analysis. If the weight factors are set 
properly, they will theoretically remove all signals that are 
coherent between the different sample locations. The goal of 




Figure 3 Multi-Channel Transversal Filter 
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such a system is to remove the correlated background noise 
from a single sensor thus enhancing signals that are 
particular to that sensor alone. If the estimated value of y k 
can be made to represent the ELF background noise, then the 
estimate error, e k , can be assumed to be due only to the 
signal present in the desired signal. The resulting output 
from the noise canceler has a much lower noise floor due to 
the removal of the correlated noise. 

2 . The Mean Square Error 

The filtering process must have a means by which an 
optimal solution can be defined. The mean-square error will 
be used to develop a defining equation, called the performance 
surface, for the optimal solution. Combining Equations 3 and 
4, a definition of estimated error in terms of the input 
vector, the desired signal, and the weighting factors is 
obtained. The expected value of the squared estimated error, 
£(e k ) 2 , assuming that e k , d k , and X k are statistically 
stationary is given as: 

E[ e{] = E[d£] +W 7 \E[X Jc X*]W-2£[d Jt x£] W (7) 

The variables x k and d, are assumed to be dependent, and in 
this case both contain the same ELF background. The equation 
for the mean square error, MSE, associated with the estimate 
is thus defined in terms of the auto-correlation function of 
the input vector and the cross -correlation function between 
the input vector and the desired signal. 
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The auto-correlation and the cross-correlation 



functions will be represented by R and P respectively. The 
elements of these two matrices are constant values if the 
input vector x k and d k are stationary random variables. A 
mathematical expression for the auto-correlation matrix, R, 
is : 



*1* • • • Vu 

*2Jfc*l* *2* ••• *2 A 



W lie • * • . 

and the cross-correlation matrix, P, is: 



( 8 ) 



P = = E[d k X lk d k X 2k . . . d k X Lk ] T (9) 

Using these two definitions in Equation 7 results in the 
following definition of the mean square error in the k ch 
estimate of y: 



MSE = J ee = E[e k ] =E[d k ] + W r RW - 2P t W (10) 

The optimum solution for the selection of the weighting 
factors is the condition where the mean square error is at a 
minimum. 

3 . Minimizing The Error 

The minimum mean square error condition can be solved 
for by taking the gradient of the MSE with respect to the 
weighting factors and setting it equal to zero. The use of 
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the minimum value of the MSE as the optimal value is based on 
the assumption that the weighting factors are only designed to 
model that portion of the input vector that is correlated with 
the desired signal. If the error is at the minimum value, 
then the filter is removing the maximum amount of correlated 
data between the input and the desired signals. The system 
being evaluated is based on the assumption that only the 
desired signal contains both noise and a signal of interest 
while the input vector is composed of a pure noise signal 
alone. Under this assumption, the only correlated data 
between the two input data samples is the noise. The desire 
is to remove as much noise as possible and thus to remove the 
maximum amount of correlated data between the input vector and 
the desired signal. It is therefore assumed that the filter 
is working at its optimum value when the weicrhting factors 
meet the criteria of minimizing the mean square error of the 
output . 

4. The Optimum Weighting Vector 

The optimum weight factor will be represented by the 
symbol W J . The two conditions that must be met in order for 
the value of J to be minimized are: 

= 0 ( 11 ) 



and 
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> 0 . 



( 12 ) 



dwjdWj 



These two equations state that when the weight vector is at an 
optimum the performance function is at a local minimum since 
the slope is zero and the curvature is upward. The gradient 
of the performance function is evaluated to be: 



V = V w J=V v E[d£] -2V w (W T P) + V w (W r RW) 
=-2P+2RW 



(13) 



Setting the gradient to zero and solving for the optimum 
weighting factors results in: 



W° = R _1 P (14) 

This equation is sometimes referred to as the Weiner-Hopf 
equation where the weight factor is referred to as the Weiner 
weight vector. The input autocorrelation matrix, R, can be 
inverted if it is "positive definite," i.e., if: 

V r RV > 0 (15) 

where V is any non-zero vector. When the matrix R is "positive 
definite" then the optimum weight factors can be solved for 
directly if the complete statistical properties of the auto- 
correlation and cross-correlation functions are known. The 
statistical properties are not known in this case so further 
development is required to obtain an estimate of the optimum 
weight factors. 
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The second condition that ensures that the error 
estimate is at a minimum is that the partial derivative of the 
gradient with respect to the weighting factors is greater than 
zero. To prove this the we take the partial derivative of 
Equation 3 with respect to W : 



a 2 j- 

dw^dwj 



= 2R 



(16) 



This result again shows that the input autocorrelation matrix 
must be "positive definite" in order to solve for the optimum 
weighting factors. 

The performance surface can now be expressed in terms 
of the minimum error possible, the weighting functions, and the 
autocorrelation of the input vector. The minimum mean square 
error possible can be found for the case when the weighting 
vector is at its optimum value. The optimum weighting factor 
is used to solve for the minimum error by substituting W° into 
the minimum error equation. 



*7*1 n = E[d z k } * [V°] T RV°-2V T V° 

= E[d k ] + [R' 1 ?] T RR' 1 P-2P T R- 1 P 
= E[dl] - P^P 
= E[dl] - P T W° 



The minimum error can now be used as the base line by which 
all other errors are found. The error at any time is the 
minimum error with the addition of uncertainty in the 
components used to derive the minimum error. Thus: 
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J = J m in + (W W°) t R (W -W°) 



( 18 ) 



The expression can be shown to equal that of Equation 10 by 
the following proof. Using the fact that (A-B) T in general 
can be expressed as A r -B : 



since all terms are scalars, the transpose is equal to the 
scalar value, thus the last two terms are equal. Substituting 
for the minimum error and the optimum weighting vector: 



This result validates Equation 18 by showing that it 
equivalent to Equation 10. An important note about this 
expression of the performance surface is that it is quadratic 
with respect to the weighting vector. 

5. Iterative Calculation of W k 

The optimum weight factors can be arrived at 
iteratively by use of a gradient search of the performance 
function. What this means is that after each estimate of the 
optimum weight vector, W K , the slope or gradient of the 
performance function, in the direction of the optimum 
weighting vector, is used to derive the next weighting vector, 
W k+1 . This method ensures that each estimation of the optimum 



J = J_ in + [W°] r R\t° +W r RW -W r RW° - [W°] T RVf 



(19) 



J = E[d k ?] -P^Pt [W°] r RW° +W r RW- 2W r RW° 

= E[d*] -P 7 )R' 1 P + P r R' 1 RR' 1 P+W r RW-2W r RR- 1 P 
= E[d*] +W T RW-2W r P 
= E[d *] +W r RW-2P r W 



( 20 ) 
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weight factors is at least as close or closer to the optimum 
value than the previous estimate. Multiplying Equation 13 on 
the left by 54 R 1 and combining the result with Equation 14, 
we obtain: 



W° ^W-Ir^V (21) 

2 

To convert this into an iterative process, the optimum 
solution is considered to be the value of the weighting vector 
for the next iteration. 

(22) 

This equation is referred to as the Newton Method algorithm 
for determining the weighting vector. 

The algorithm takes only one iteration to arrive 
theoretically at the optimum solution for the values of the 
weighting vector. This algorithm would be the ideal method to 
use if exact solutions to the gradient term and the inverse 
correlation of the input vector were known exactly. This is 
not generally the case and an estimate for the gradient must 
be used. To compensate for the lack of knowledge of the 
gradient at iteration k, a constant p, will be used instead of 
the factor Vt in Equation 22. The requirements on the range of 
values that p can take on are: 
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0 < |i < 1 



( 23 ) 



The effect of p on the solution at each iteration is that it 
affects the rate of convergence of the algorithm. A 
discussion of the attributes and further limitations on the 
value of p will be covered later. A general form of the 
Newton's method of gradient search is therefore: 

= 1 ( 24 ) 

This expression is ideal under the following three conditions: 

1 . \1 = V2 

2. Exact knowledge of the gradient vector, V k . 

3. Exact knowledge of the (unchanging) matrix R 1 . 

These conditions, unfortunately, are not normally 
attainable. The value of p is selected between 0 and V* so 
that the algorithm is overdamped and can accommodate 
fluctuations in the matrix R. The compensation for this 
selection of p is that the algorithm takes longer to obtain 
its optimum solution. A second modification to the value of 
p is to base its value on the average eigenvalue, A av , so that 
p is replaced by pA, av . 

The second condition pertaining to the exact solution 
to the gradient of the error surface is, in general, never 
attainable. The gradient must therefore be estimated in some 
manner. The estimate of the gradient that will be used is 
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based on the square of the present value of the error signal 
and not the true ensemble average as it has been up to now. 




(25) 



From the definition of the error signal: 



d *k _ 



(26) 



aw* aw* 



which leads to the following estimate of the gradient: 



With these two relaxations from the three ideal conditions, 
the algorithm is now in a form known as the Newton/LMS 
algorithm [Ref 7]: 



This is an optimum iterative algorithm that is only complete 
if full knowledge of the input correlation function is known. 
This is generally not the case and thus a final modification 
must be done to produce a working estimate to this algorithm. 
The final modification is an iterative estimate of the inverse 
of the input correlation matrix, R" 1 , and is called the 
Sequential Regression (SER) algorithm. 



^ k - ZtiXk 



(27) 



W * +1 = W^/i^R- 1 ** 



(28) 
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IV. THE SEQUENTIAL REGRESSION (SER) ALGORITHM 



A . DEVELOPMENT 

The Sequential Regression (SER) algorithm uses the 
LMS/Newton algorithm as its basis in conjunction with a method 
for iteratively estimating the value of the inverse input 
auto-correlation matrix to produce the weighting vector for 
the filter coefficients. The development of the algorithm 
comes from Adaptive Signal Processing by Widrow and 
Stearns [Ref 7], The core of the algorithm is the use of an 
estimate for the inverse of the input auto-correlation matrix 
which reduces computational load with minimal estimation 
error . 

The input auto-correlation matrix is 

R =E[X k xl] (29) 

where the subscript k goes over the entire length of the 
random process X } . R will be estimated iteratively using only 
a finite or truncated sample of the random process X,. The 
estimate of the input auto-correlation matrix is given by: 

= -A- t X jXj. (30) 

k k+1 j-o 1 

This estimate is unbiased when the input variable, X lf is a 
stationary random process. When X l is not stationary, it can 
be shown that the estimate is not very accurate. The estimate 
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of R.. at each value of k is equally influenced by the 
previous, k-1, estimates of the input auto-correlation matrix . 
In order that the most recent estimates of R k have more 
influence over the new estimate, a "fading" memory term, a, is 
introduced which reduces the significance of the past 
estimates. The fading memory expression of the input auto- 
correlation matrix is denoted by Q v where: 

Q* = (31) 

The value of a is chosen, as a rule of thumb, such that the 
half life of the exponential function is equal to the length 
of stationarity of X, [Ref 7] . This rule of thumb leads to the 
following statements about the value of a. 



The value of the fading memory factor over k iterations is: 



and thus the estimate of the input auto-correlation matrix 
becomes : 



g K 2 (-1/length of stationarity of X) 



(32) 



0 < a < l 



(33) 




(34) 




k 



( 35 ) 
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This estimate is exact for the condition where X K is constant. 
If the input variable, X ; , is stationary, a is one, and if the 
limit of Equation 35 is taken as a approaches one, the result 
is in agreement with Equation 30. 

The derivation of the SER algorithm begins with the 
estimate of the input auto-correlation matrix. Using this 
estimate in the equation for the optimum weighting factors 
results in 



Mi = **• ' 36 > 

The definition for the estimate of the cross-correlation 
vector, P k , is, similarly to Equation 35, 



= 



1-tt 

l-a ^ 1 



a 



k-1 



d,Xj. 



( 37 ) 



Using Equations 35 and 37 in 36 gives: 

Q > cWjc = (38) 

1=0 



The SER algorithm attempts to estimate the next value of 
the weighting vector, W K+1 , based on the present values of R,, 
and P v . To do this. Equation 38 becomes: 
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( 39 ) 



9A.. - I> ,k -"' 1 d 1 X 1 +d k X k 

i=G 

= a Qk-l^ k + 

= (Q k -Tp2>* k *dp k . 

The last line of Equation 39 is the result of the following 
relationship : 

Qk = «Qk-i + *k*l (40) 

The value of d k can be replaced by recalling Equation 5 and 6 
which state that 



<4 = *k + y k = 

thus Equation 39 becomes 



( 41 ) 



«A., = * (**♦ *!"*>** 

( 42 ) 



= 0) Fk + e k*k‘ 

Multiplying on the left by Q*' 1 , gives an iterative method for 
calculating the weighting factor vector 



w*.i - * 0*t*x*- <«) 

The expression for W k+1 in Equation 43 is similar to the 
LMS/Newton method derived earlier. With this in mind, an 
approximation for the LMS/Newton algorithm can be made using 
the above derivation. The SER algorithm used as an 
approximation of the LMS/Newton method is given as: 
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(44) 



"*♦! = W ic + 



2 ^ave( 1_aiC+1 ) 

1-a 



Qi le Jc X /c‘ 



The algorithm requires a method of iteratively solving for the 
value of Qv 1 in order to be complete. 

The method to iteratively solve for Q,/ 1 begins by pre- 
multiplying Equation 40 by Q k ' x , post-multiplying by Q^' 1 , and 
then by X k , to obtain: 



Q'kQt&k- iXjc = Q^aQjc-iQi-iXjc + 
which reduces to, 

Qi- iX* = *Q*x k +Q~kZi?ZQ£iX k 

= Q~ k 1 X k (a ^xlQ~ k 1 . 1 X k ) 



(45) 



(46) 



The quantity in parentheses is a scalar and can be divided out 
while the quantity X k T Q k _ 1 ' 1 is also multiplied on the right to 
give 



a +xj&£ 1 x Jt 



QiXjtXSsii • 



(47) 



The first lines of Equation 46 and 47 can be combined and re- 
arranged to form. 



Q'* 1 




(Q^-iXj (o^ix^r 

*+Xk(Qi-i* k ) 



(48) 



An iterative method of calculating Q k _1 is now available and is 
only based on its previous value, Q^' 1 , and the present value 
of the input, X k . The value in parentheses is used three 
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times, it will be calculated separately and represented by 
S k . 

The SER algorithm is now completely derived as an 
iterative method of approximating the LMS/Newton method of 
adaptive filtering. The following summary shows the process 
by which each new value of the weighting factor vector is 
calculated . 



g 2 (-l/lengthof etationaiity of X) 
e k ~ 

5 = Qi-xX k 

Y = a+x£s 





1 



- w* + 



2 \i\ ave {i-a k * 1 ) 



l-a 



Q* z k*-k 



o < ^ < *r~ — ; 



( 49 ) 



The values that must be carried forward from one iteration to 
the next are Q^' 1 , and W k . These represent the previous value 
for the estimate of the inverse input auto-correlation matrix 
and the present estimate for the weighting factor vector. 

The output of the algorithm is the error, e k , between the 
present value of the desired signal and estimate of the 
desired signal: 
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( 50 ) 



OUTPUT = e k = d k -W k Xl 
This value is saved after each iteration in a vector 
representing the uncorrelated portion of the desired signal 
and noise vector of the primary channel with the noise only 
vectors of the reference channels at the time index k, 

X k . [Ref 7] 

B. INITIALIZATION 

The SER algorithm assumes that the value of the inverse 
input auto-correlation matrix and the weight factors vector is 
known from the previous step. At some point in the past, an 
estimate for these terms must have been made so that a 
starting point can be established. The required accuracy of 
the estimates depends on the length of time that the algorithm 
will be operating and the earliest time that accurate data is 
required . 

The driving force behind the accuracy of the initial 
estimate is the speed with which the algorithm can converge to 
an optimal solution. The speed of convergence is dependent on 
two elements of the algorithm. The first element is the 
initial value of the inverse input auto-correlation matrix, 
Q k _1 . The significance of this initial element is that the 

further the initial estimate is from the true value, the 
further the algorithm has to adapt in order to approximate the 
optimum solution. The second element is the step size, 2pA, ave , 
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which specifies how far the algorithm can progress towards the 
optimum value in one iteration. The step size will be 
discussed in the next section but for now it shall be assumed 
to be a constant in this discussion of the initial conditions. 

The use of a large constant multiplied by the identity 
matrix as the initial value of the Q k ' a has been sufficiently 
argued by Lee [Ref. 8] . 

Qo 1 = g 0 i (51) 

As discussed by Lee, the use of a large constant for q 0 is the 
proper selection for the case where the random process is 
stationary. Widrow and Stearns argue that this same choice 
for the initial value of q 0 , in a non-stationary condition, is 
also appropriate [Ref 7], The point is made that if a priori 
knowledge of the random process allows for a better estimate 
of the value of Qo" 1 , then that information should be used. 
The method of initialization for this study is to use Equation 
51 with q 0 equivalent to 100. 

The initial values used in the weighting factor vector are 
assumed to be zero. This assumption is made since knowledge 
of the process before time zero is unknown. The use of a 32 
tap delay filter also allows for an estimate of the initial 
full weight factor vector but this is done so at the expense 
of the first 32 output data points. This is done by expanding 
the number of elements in the weight factor vector as each of 
the first 32 points are read into the algorithm. The 
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algorithm is applied only to the data that is available; 
therefore, the first time that each element of Vl r is non-zero 
is at time index 31. This corresponds to the time when the 
first 32 data points are input into the algorithm. The result 
of this initialization is that the first real data point to be 
considered an output of the SER method used here is at k = 31. 
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V. STABILITY AND CONVERGENCE 



A. THE SINGLE WEIGHT CASE 

The SER algorithm contains values that can be set to 
affect the rate at which the algorithm achieves the optimum 
solution. A simple case where only one weight factor is 
concerned will be used to demonstrate the relationships 
involved in the selection of the convergence rate factors. 
The single weight iterative process using a gradient search 
method like the SER algorithm can be represented by: 

W k+1 = + ( 52 ) 

The expression for the single weight gradient at time 
index k is given by: 




d[\ (w- w°) 2 ] 



= 2\ (w k -w°) 



( 53 ) 



The value of the input auto-correlation matrix, r 00 for the 
single weight case, is replaced by its equivalent eigenvalue, 
A,=E[x£] . Combining Equations 52 and 53, an iterative equation 
for the weight factors is found: 
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***♦1 = w k -2\i\{w k -w c ') (54) 

for the single weight case. The terms in Equation 54 can be 



rearranged to combine like terms. 



w k =w°+ (1-2|aA) k {w 0 -w°) 



(55) 



This expression contains the geometric ratio of successive 
terms which is also used to define the stability requirements. 

Stability is guaranteed if the absolute value of the 
geometric ratio is less than one, i.e., 



The limits placed on the value of the u can be expressed in 
terms of the eigenvalue of the input auto-correlation matrix 
X by rearranging Equation 56 



The optimum solution can be seen as the solution of Equation 
55 if the limit of w k as k goes to infinity. This is because 
the geometric ratio approaches zero as k approaches infinity. 
The rate at which the algorithm approaches the optimal value 
is dependent on the value of the geometric ratio. 

Figure 4 shows the effect of varying values of the 
geometric ratio on the rate of convergence. If the value of 
the geometric ratio can be maintained between zero and one, 
then the algorithm will always approach the optimal solution 
from one side and can be considered as overdamped. For a 



|l-2jiX|<l 



(56) 




(57) 
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geometric ratio of zero, the algorithm reaches the optimum 
value in one iteration from one side and can be considered as 
critically damped. If the magnitude of the geometric ratio 
falls between zero and negative one then the algorithm is 




Convergence: overdamped 0<g<l; critically damped g=0; 
underdamped -l<g<0 



B. MULTI -WEIGHT CONDITION 

The single-weight condition discussed in the previous 
section can be expanded to develop an understanding of the 
multiple weight factor case. The same conditions that applied 
to the single weight case are true in the multiple weight 
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case. The difference is that each weighting factor has an 
associated eigenvalue and thus the condition of stability is 
now dependent on a set of simultaneous equations that must 
meet the stability condition. 



W Jc =W°+ 



(l-2|iA 0 ) 

( 1 - 2 ^) 

(l-2gA 2 ) 



(W 0 -W°) 



(1 2\iX L _ 1 ) 



( 58 ) 



The subscript on A indicates the index value corresponding to 
a specific weight value. There are a total of L weighting 
factors and thus there are L individual eigenvalues, A. 

The value of g in each of the simultaneous equations is 
held constant. Since g is the same for all equations, the 
eigenvalue that has the largest value will impose the most 
restrictive conditions on the selection of a value for p. The 
new range of values that p can take on is: 

^->^>0 ( 59 ) 

^max 

The condition stated in Equation 59 gives the multi-weight 
stability requirements that must be met to ensure that the 
algorithm will converge to the optimal solution. 

The rate at which the algorithm will converge was shown in 
Figure 4 to be dependent on the value of the geometric ratio. 
Now that the value of p is a function of the maximum 
eigenvalue of the input auto-correlation matrix, R, it can be 
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shown how this affects the rate of convergence. By holding 
the value of p constant for each weighting factor, an 
undesirably slow convergence rate can occur for widely varying 
values of 1. As shown previously, the algorithm converges in 
one iteration when p is equal to 1/2A, for the single weight 
case. Using that same criteria, the value of the multi-weight 
constant p is l/2X max . The use of this criterion ensures that 
the algorithm falls into the overdamped category since p is 
always less than or equal to 1/2^ for 1 = [ 0 , 1, 2, ..., L] . 
The weight factor, w : , with the smallest value of k will have 
the slowest convergence rate since its geometric ratio will be 
farthest from zero. This is true because the value of 2\ik 
will be at its minimum and will cause the geometric ratio to 
approach one. When this condition holds, the algorithm may 
take excessively long periods of time to converge to an 
optimal solution. 

C. NON- STATIONARITY EFFECTS 

The discussion up to this point in dealing with stability 
and convergence assumed that the process under study was 
stationary. Under stationary conditions, the values of the 
eigenvalues of R are assumed to be constant. A second key 
assumption to the development is that the values of the 
eigenvalues are known. The lack of knowledge of the input 
auto-correlation matrix, R, requires that some estimation be 
made for the value of p. The value of p must be adjusted to 
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ensure convergence of the algorithm in an environment where 
the eigenvalues of the auto-correlation matrix are varying at 
some unknown rate from iteration to iteration. 

A method for compensating for the variations in the range 
of values that p can assume is to use an averaged eigenvalue, 
A. ave , for the selection of the constant p for each iteration. 
The value of p is now a constant for a single iteration but is 
changed from iteration to iteration based on the average 
eigenvalue. Stability of the algorithm is guaranteed if the 
value of the product, pA. ave , is between zero and one. It is 
important to note that the expression pA ave is now used only to 
replace the constant p used previously. The idea is to scale 
the average eigenvalue for that iteration by a constant, p, 
and use it as the factor for controlling convergence of the 
algorithm. 

The geometric ratio which controls the convergence of the 
algorithm has now been altered to include a term to compensate 
for the non-stationarity of the input random process, X. The 
new expression for the geometric ratio is 

|1 - 2pA ave R- 1 | < 1 (60) 

The inverse input auto-correlation matrix, R 1 , is included in 
the development of the SER algorithm. The convergence 
controlling term is thus pX ave , and must be less than the 
maximum eigenvalue of R to ensure stability and convergence. 
A second condition that can be imposed on the product pA, ave is 
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that its value must always be between zero and one. Using the 
fact that the spread between maximum and minimum eigenvalues 
can be rather large, the product of pA ave must be much less 
than one . 

The ELF background noise processes have shown a tendency 
to have a wide range between eigenvalues and to slowly change 
over time. Because of these two conditions, the value of the 
pA. ave must be adjusted over time to account for these 
variations. The chosen method is to estimate the average 
eigenvalue by taking the average of the square of the input 
vector, X k . This value is multiplied by the 1/100 of the 
inverse of the maximum of the squared input vector. 



\*^ave 



1 

100 [max (Xj) ] 2 




( 61 ) 



The use of this estimate for convergence control in the 
algorithm ensures that the requirements for stability are met 
for each iteration. The 1/100 term ensures that the maximum 
value of p is not exceeded. It was derived empirically with 
the use of an actual set of ELF background measurements and 
may require adjustment in a different environment. 

The choice of using the present input to the algorithm to 
get an estimate of the average eigenvalue is used for 
computational ease with minimal added error. The use of the 
inverse input auto-correlation matrix should provide a more 
accurate measure of the eigenvalues; That, however, requires 
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the inversion of the matrix. The prevention of matrix 
inversion is the exact reason the SER algorithm is developed, 
so to obtain the eigenvalues in this way defeats the reason 
for the algorithm’s design, and is therefore counter- 
productive. The method has proven adequate and reasonably 
accurate in empirical work up to this point, but does leave 
room for improvement and further study. 
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VI . IMPLEMENTATION 



A. A PHYSICAL DESCRIPTION 

The Applied Physics Lab has two identical sensor platforms 
submerged in shallow water. Each platform has electrode pairs 
for detecting the geoelectric ELF background noise. The 
system has been in operation for over two years collecting 
samples of the background noise environment. 

Ideally, the two platforms are parallel when the platforms 
are placed on the bottom of the bay. Figure 5 shows the ideal 
alignment of the platforms. The two platforms are labeled as 
the east, E, and west, W, platforms. Each platform has a set 
of orthogonal electrode pairs, labeled X and Y. These 
electrode pairs are orientated to receive the two orthogonal 
source vectors, E x and E y , in the horizontal plane. The 
strongest geoelectric field strengths are in the horizontal 
plane due to the shift in polarization as a wave crosses the 
water-air boundary from vertical to horizontal. 

If the two platforms are perfectly aligned and the noise 
recorded on each are identical, the noise in a specific 
direction, X or Y, can be removed by subtraction. To cancel 
the x-direction noise simply subtract XW from XE. Because the 
platforms are not perfectly parallel and the recorded noises 
are not identical though highly correlated, this method is not 
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Figure 5 Illustration of Sensor Set-Up 

ideal. The actual platforms can be misaligned up to 10 
degrees from perfectly parallel. This is due to physical 
constraints such as bottom contours and mounting technique. 
The result of such a misalignment is that both the XW and YW 
sensors of the west platform have components that are coherent 
with XE of the east platform. If the electrode pairs in 
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Figure 5 are thought of as vectors, then it is easily seen how 
misalignment can cause coherency between XE, XW, and YW. 

A set of real data is traced through the signal processing 
system designed to accompany this sensor system. A calibrated 
source signal is embedded in the signal XE. The path that the 
received signal follows can be seen in Figure 6. The output 
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Figure 6 The ELF Signal Processing Block Diagram 
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of each sensor is sampled at 128 Hz. The data is bandpass 
filtered and resampled down to 32 Hz for processing. The data 
is then passed through the SER algorithm where the coherent 
portions of the west platform are removed from XE . The output 
of the SER portion of the processor is demodulated and passed 
through a likelihood ratio detector where a decision can be 
made as to when the embedded signal is present. 

B. INPUTS TO THE SER ALGORITHM 

The SER algorithm is used to reduce the ELF background 
noise detected in the primary submerged electric field sensor 
by cancellation with the output of the reference parallel 
electric field sensor. These parallel reference sensors are 
also submerged and are, ideally, assumed to be located within 
the same ELF background noise as the primary channel. The 
primary channel y represents sensor XE and the reference 
channels, xl and x2, represent XW and YW, respectively. The 
samples are the actual output of the sensor system used by 
APL. The samples are taken at the same time without any known 
signals present. The sampling rate of the data input for the 
SER algorithm is 32 Hz and the total length of data is 113760 
points. These data samples correspond to almost an hour of 
background noise measurements . 

The filter is a 32 point transversal filter that is 
applied to each channel simultaneously. The length of the 
transversal filter means that the previous 31 samples of the 
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sensor are used with the present sample value for each 
iteration of the algorithm. The goal of the algorithm is to 
take the output of the two reference channels, xl and x2, and 
to generate an estimate of the noise in the primary channel y. 
Figure 7 shows what the primary channel's data looks like in 
real time. The length of stationarity , used to determine the 
forgetting factor, a, is assumed to be about 300 data points 
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Figure 7 The Primary Noise (Y) and Calibrated Source Signal 

long. This value is used based on empirical analysis of the 
real data. 

The effectiveness of the algorithm is analyzed by running 
a set of scenarios with a known signal inserting in the data 
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samples. The scenarios consist of manually inserting a 32 Hz 
signal on a 5 Hz carrier that is 8192 data points long. 
Figure 7 shows the components of the primary channel signal in 
a real time representation. The primary channel is processed 
by the SER algorithm while varying the strength and location 
of signals embedded in the background noise y. 

C. SIGNAL STRENGTH RESULTS 

The algorithm is first run with the location of the signal 
within the noise held constant. The strength of the signal is 
varied by dividing each data point in the signal vector by a 
constant. A series of four runs are conducted, each 
succeeding run has half the signal strength of the preceding 
run. The initial run uses the signal shown in Figure 7 with 
the amplitude divided by 500. The signal is indistinguishable 
prior to processing as seen in Figure 8. 

The signals for each run are centered at time 19.2, 32 and 
44.8 minutes into the noise sample y. There is no 
significance to the positions chosen except that the signal 
data streams were not allowed to overlap. The objective of 
this first set of runs is to demonstrate how the SER algorithm 
enhances the ability to distinguish small signals buried in 
the noise. The data is first processed in the SER algorithm 
to remove the coherency between the reference channels and the 
primary channel. The output is then demodulated leaving only 
the in-phase and quadrature components of the non-coherent 
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portion of the primary channel. An identical reference set of 
data is demodulated without passing it through the SER 
algorithm so that a visual and quantitative analysis of the 
effects of the SER algorithm can be made. The output of the 
demodulator is shown in Figures 9 and 10. Note that the 
envelope of the signal is visible to nearly the strength of 
l/2000 th , or the third run made. 

The output of the demodulator is passed as input into a 
matched filter detector based on the demodulated signal. The 
output of this detector is shown in Figures 11 and 12 with the 
SER and non-SER output for a given signal strength shown 
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Figure 8 Input Primary Channel Noise And Signal 
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together. Numerically, the peak value out of the detector is 
three to five times greater when the SER algorithm is used as 
compared to when it is not used. 

A final numerical analysis is performed by using a primary 
channel without a signal installed to get a measure of the 
noise present in the output. The standard deviation of this 
output of the detector when this noise only vector is used as 
input is used as an approximation to the noise power. 
Dividing the outputs of the detector in Figures 11 and 12 by 
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Figure 9 Demodulator Output of 
(a) Signal/500 (b) Signal/1000 
(d) Signal/4000 
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Figure 10 Demodulator Output of Quadrature Components; 
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Figure 11 In-Phase Detector Outputs: (a) Signal/500 

(b) Signal/1000 (c) Signal/2000 (d) Signal/4000 
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Figure 13 SNR vs Relative Signal Strength 



D. DETECTOR THEORY 

The process of adaptive noise cancellation is performed 
for the purpose of increasing the signal-to-noise ratio. A 
section on the narrowband detection schemes that have been 
considered in conjunction with the noise canceler are 
presented since they are integral to the analysis here. Two 
likelihood ratio detectors are considered for the detection of 
narrowband signals . 



54 



The first detector is designed with the assumption that 
not only is the signal known, but its peak within a given 
interval is constant relative to the start of that interval. 
This last assumption is to assure that the demodulated in- 
phase and quadrature signals that are used in the likelihood 
ratio detector match those in the noise canceled and 
demodulated output of the sensor system described previously. 
In the second detector, the latter assumption was dropped and 
the design followed that of an optimal envelope detector. In 
both cases it was assumed that the residual noise field was 
gaussian. This assumption, while not strictly true, seems to 
yield quite reasonable results . [Ref . 9] 

The signal is generated by a calibrated source with known 
qualities. The modelling equation for the narrowband signal 
is given by: 

s(t) = A { t ) COS [ 2 7i f c t + 0 ( t ) ] (62) 

For the first detector, the phase term, 0 , is assumed constant 
within any non-overlapping interval 0 to T. The interval of 
0 to T is equivalent to the length of the signal. The purpose 
of this assumption is so that a Matched Filter can be created 
for the in-phase and quadrature components by using the 
demodulated in-phase and quadrature component of the signal. 
Figure 14 shows the demodulated signals used for the detector. 

The first detector is a multi-channel Matched Filter shown 
in Figure 15. The Matched Filters, h i and h ql are determined 
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Figure 14 Demodulated Signal for Match Filtering 

based on the requirement that the output signal-to-noise ratio 

is maximized at time T. The corresponding equations are 



/ 



R ii ( t-x ) R iq (t-x) 
R ql (t-X) R qq ( t-x) 



h d (x) 
h g (t) 



dx = 



Si {T-x) 
Sg{T~X) 



(63) 



where the quantities R denote the various auto and cross 
correlations. This detector is the one used for all 
processing performed in this thesis. The input to the 
detector is the noise cancelled output of the SER portion of 
the processor after demodulation . [Ref . 10] 

The second detector is an optimal nonlinear envelope 
detector. The description of its principles and design are 
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located in the Appendix. The detector was not used by the 
author in this work relating to the analysis of the SER 
algorithm, but it is presented as complementary work relating 
to the overall project. 

E. SIGNAL LOCATION 

The signal strength analysis showed that the signal 
divided by 2000 is easily discernable from the background at 
the output of the detector. The effect of varying the 
location of the signal will allow analysis to get a better 
feel for the expected SNR based on a larger sample size for 
averaging than previously used. A second benefit of this 
analysis is to see the effect of the initialization process on 
the SNR. The signal will be shifted by a quarter of its 
length for each run, or 2048 points. The SNR is plotted 
versus the location of the center of the signal in Figure 16. 

The data is run with and without the SER algorithm to get 
a good feel for the enhancement achieved with the SER 
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In — Phase 





Figure 16 SNR VS Location of Signal In Noise Vector (Y) 

algorithm. Three separate locations of the signal will be 
used for each run to reduce calculation time. The use of 
multiple signals causes no effect on the SNR calculation as 
long as the signals are not allowed to overlap. 
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VII. RESULTS AND CONCLUSION 



A. ANALYSIS OF RESULTS 

The use of the SER algorithm enhances the ability to 
detect narrowband low power signals buried in the background 
noise. Figure 17 shows the spectrum of the noise sample from 
sensor XE before and after the use of the SER algorithm. The 
frequency range is from one to ten hertz, which is in the 
range of concern for the JHU project. The spike seen at four 
hertz (note the x-axis is a logarithmic scale) is caused by 
the processing equipment. This is a man-made noise source and 
will be filtered out in the final project design. The figure 
illustrates the 10-15 dB of noise cancellation possible by 
using this method of adaptive noise cancellation. 

The location of the signal in comparison to the start of 
the adaptive process has been shown to have very little effect 
on the signal enhancement. The system implementation will 
most likely cause the algorithm to only require initialization 
once a week and operate constantly otherwise. The 
initialization process takes about 100 iterations or 1.67 
minutes before settling out but is heavily data dependent. 
Even while settling, the performance of the algorithm is 
adequate to receive most signals. For the purpose of this 
project, initialization should be of no concern. 
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Figure 17 Illustration of SER Algorithm Cancellation 
Capability, Power Spectral Density of 1-10 Hz 

B. CONCLUSIONS AND RECOMMENDATIONS 



It was hoped at the beginning of the research that 
information could be gained on the characteristics of the 
noise environment . The stationarity of the background noise 
environment was found to be at least nine seconds. This is 
believed to be only a characteristic of the specific data set 
that is used. The length of stationarity affects the value of 
the forgetting factor, o, and is implied by the use of the SER 
algorithm. A better feel for a good average value to use for 
the length of stationarity will result from the use of more 
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data sets. This analysis needs to be conducted and could be 
a source of further work in this area. 

A second area of work that could develop from this 
research is a study of the cancellation possible from a 
vertical above-surface sensor and a submerged horizontal 
sensor. The purpose of this would be to differentiate between 
surface and submerged activity. This may prove to be 
impossible but it would be interesting to see the amount of 
correlation possible between the two locations. This type of 
arrangement could also lead to a wider distance between 
sensors which could also increase the knowledge of the 
background environment . 

The purpose of the research presented in this thesis is to 
develop a method to adaptively cancel the ELF noise present in 
a submerged sensor array. The SER algorithm appears to be an 
adequate method to accomplish this goal of noise cancellation . 
The cancellation ratio in the frequency range of interest is, 
on average, 12.5 dB. This cancellation ratio along with the 
use of optimal detectors will assist the system in detecting 
narrowband signals buried deep in the background noise. 
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APPENDIX 



The second detector used in conjunction with this project 
of measuring the geoelectric background noise is an optimal 
nonlinear envelope detector. The symbolism used is slightly 
different but efforts are made to define all relevant terms. 
The in-phase and quadrature data (after noise cancellation) 
are denoted by and r^, respectively. Similarly the in-phase 
and quadrature noises are represented by nj) and n£. Assuming 
a joint Gaussian distribution we have the following two 
hypothesis : 

H 0 : ( r 1 , r q ) = ( n 1 , n q ) 

H. : ( r 1 , rM = ( , n q ) + ( s 1 , s q ) 

where H 0 denotes the null hypothesis (noise only) and H x 
denotes the alternative (signal and noise present). The 
vector notation means 



r = [r 1 , r 3 ] = [ (r* , rf) , . . . , (r£,rg)] 
where N denotes the length of the interval for which the 
decision is being made. It is then clear that the probability 
density function for the received data is as follows: 



P (r|«i) 




( r-s ) t R~ 1 ( r-8 ) 

— nn 



where E rr , denotes the covariance matrix of the noise vector. 
In this particular case, it is assumed that the in-phase and 
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the quadrature components have been decorrelated and whitened 
with the same variance o so 



p<*I«l) = exp f^rr£ (^*-s*) 2+ 2 

\ 2na 2 ) { 2 a 2 [ti pi 

Introducing polar variables ( , <p\ ) and ( rf , <f% ) it is 

found that 



N 



II r * 




N N N 


P = ^^5 exp< 

27 UJ 2 


1 

2o 2 


£ ** + £ S l~ 2 £ r * S *COS (<}>*) ' 

Jc=l k = 1 Jc=l J 



where r 2 = r * 2 + r % 2 and tan ( <|>, ) = . Integrating 

the above over the angular variables will give 



N 




where T 0 is the modified bessel function of order zero. The 
vector r on the left hand side now only refers to the various 
envelopes, r = [ r lt . . ., r N ] . The same is true for the signal 
envelope. The likelihood ratio detector is found to be: 



P(r\H 0 ) 



=exp 





The log-likelihood ratio is then: 



N N 

logA=- — -£ s* + £ log 

2o 2 £=i Fi 
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The first term is a constant, dependent only on the known 
signal, and can be dropped. Finally the test statistic or 
detector output is [Ref 10]: 

logX-glog |l 0 (f^i|. 

The implementation and analysis of this detector scheme is 
left as a possible follow-on thesis topic. 
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